/*
 * CanonicalReactionSystem.cpp
 *
 *  Created on: 7 Sep 2011
 *      Author: Allan
 */

#include "CanonicalReactionSystem.h"

// C++ includes
#include <cmath>
#include <iomanip>
#include <set>

// Boost includes
#include <boost/foreach.hpp>

// Eigen includes
#include <Eigen/LU>

// GeoReact includes
#include "PrimarySelector.h"
#include "Utils.h"
#include "Utilities/Algorithms.h"
#include "Utilities/Assert.h"

CanonicalReactionSystem::CanonicalReactionSystem() 
{}

CanonicalReactionSystem::CanonicalReactionSystem(const ReactionSystem& reactions) :
ReactionSystem(reactions)
{
	PrimarySelector selector(reactions);
	
	SetPrimarySpecies(selector.SelectByFrequency());
}

CanonicalReactionSystem::CanonicalReactionSystem(const ReactionSystem& reactions, const vector<string>& primarySpecies) :
ReactionSystem(reactions)
{
	SetPrimarySpecies(primarySpecies);
}

CanonicalReactionSystem::~CanonicalReactionSystem() 
{}

void CanonicalReactionSystem::SetPrimarySpecies(const vector<string>& primarySpecies)
{
	// Set the primary species
	jSpecies = primarySpecies;
	
	// Set the secondary species
	eSpecies = Difference(GetSpecies(), jSpecies);
	
	Assert(eSpecies.size() == GetNumReactions(), 
			"The provided set of primary species yields an "
			"incorrect number of secondary species.");
	
	// Assemble the stoichiometric matrices (uej) and (uee)
	uej = AssembleStoichiometricMatrix(jSpecies);
	uee = AssembleStoichiometricMatrix(eSpecies);
	
	// Compute (inv_uee), which is the inverse of the stoichiometric matrix (uee)
	FullPivLU<MatrixXd> lu(uee);
	
	inv_uee = lu.inverse();
	
	// Compute the canonical stoichiometric matrix of the reactions
	vej = - inv_uee * uej;
}

const vector<string>& CanonicalReactionSystem::GetPrimarySpecies() const
{
	return jSpecies;
}
	
const vector<string>& CanonicalReactionSystem::GetSecondarySpecies() const
{
	return eSpecies;
}

const MatrixXd& CanonicalReactionSystem::GetCanonicalStoichiometricMatrix() const
{
	return vej;
}

const VectorXd CanonicalReactionSystem::CanonicalEquilibriumConstants(double T, double P) const
{
	VectorXd lnKe = EquilibriumConstants(T, P).array().log();
	
	return (-inv_uee * lnKe).array().exp();
}

const VectorXd CanonicalReactionSystem::CanonicalReactionQuotients(const VectorXd& a) const
{
	VectorXd lnQe = ReactionQuotients(a).array().log();
	
	return (-inv_uee * lnQe).array().exp();
}

const VectorXd CanonicalReactionSystem::CanonicalKineticRates(double T, double P, const VectorXd& n, const VectorXd& a) const
{
	return uee.transpose() * KineticRates(T, P, n, a);
}
